clear all
clc
rng(1342,'twister') 
s=rng

%% Data Prelims %%

%Load data
filename = 'D:\Projects\HeterogeneousExternalities\Data\SimulationMatrix.xls';
load = xlsread(filename,'A:L');

%Number of counties
size=380;

%Variable inputs
pd=load(1:size,1);
vmt=load(1:size,2);
eta=load(1:size,3);
beta=load(1:size,4);
speed=load(1:size,5);
v0=load(1:size,6);
VOT=load(1:size,7);
mpg=load(1:size,8);
fuel=load(1:size,9);
urban=load(1:size,10);
statetax=load(1:size,11);

%Scale problem to help solver
v0=v0/1000;
vmt=vmt/10000;

e=eta*(-1);

%Replace positive betas with zero
beta=min(0,beta);

%Problem size
n=length(e);
z=zeros(n,1);


%State tax converted to per mile tax
statetax=statetax/100;
ts=statetax./mpg;

%Fuel costs per mile
F=fuel./mpg;

%Scale vehicle counts to vmt
theta=vmt./v0;
tv0=theta.*v0;

%Find alpha using average values
alpha=speed-(beta.*v0);

%Find A using average values
A=(tv0)./(F.^(-e));
%% Integration %%

%Maximization variable
v=sym('v',[n 1]);

%Integrate over demand function
for i=1:n
    counter=i
    B(i,1)=int((A(i)/(v(i)*theta(i)))^(1/e(i)),v(i), [.1 v(i)]);
end

%% Maximiation %%

%Starting values
starter=v0-rand();
%Clear out prior optimization if necessary
v=sym('v',[n 1]);
tv=theta.*v;
P=(A./tv).^(1./(e));
t=P-F;
C=tv.*((F-ts)+pd+(VOT./(alpha+(beta.*v))));
y=C-B;
f=sum(y);
objfun=matlabFunction(f,'vars',{v});

%Revenue from "optimal" tax
opt_rev=(t+ts).*tv;
tot_opt=sum(opt_rev);

%Revenue from current taxes
s_rev=ts.*tv0;
tot_s=sum(s_rev);

%Revnue constraint
srev_con=tot_s-tot_opt;
sceq= [srev_con];
%Nonlinear inequality constraints (v greater than 0)
sc = [-v];

sconfun = matlabFunction(sc,sceq,'vars',{v},'outputs',{'sc','sceq'});

%List of maximum tax constraints
tlist=[.01;.1;.2;.3;.4;.5;.75;1;1.5;2.25;3;4;6;10;14;20;inf];

welfare=zeros(length(tlist),1);

for i=1:length(tlist)
    counter=i
   
%Adjust max tax constraint
max_tax=tlist(i,1);
tax=max_tax./mpg;

%Lower and upper bound for v---lb is v at max tax, ub is v at tax=0
lb=(A.*((F+tax).^(-e)))./theta;
ub=(A.*((F-ts).^(-e)))./theta;

%Global solver
rng(s) 
problem = createOptimProblem('fmincon',...
    'objective',@(v)objfun(v(:)),'lb',[lb],'ub',[ub], 'nonlcon',sconfun,...
    'x0',starter,'options',...
    optimoptions(@fmincon,'Algorithm','sqp','Display','off','MaxFunEvals',1500000,'MaxIter',10000));
gs = GlobalSearch('Display', 'iter');
[v eflag output] = run(gs, problem);
%Replace starting values with latest optimization run
starter=v;
tv=theta.*v;

P=(A./tv).^(1./(e));
t=P-F;
t_pg=t.*mpg;
opt_rev=(t+ts).*tv;
tot_opt=sum(opt_rev)

w_no=-objfun(v0);
w_opt=-objfun(v);
w_gain_opt=w_opt-w_no

%Convert from daily to annual welfare gains
answer=w_gain_opt*365;
vmt_opt=tv;
fulltax=t_pg+statetax;
welfare(i,1)=answer(1,1);

dir='D:\Projects\HeterogeneousExternalities\Output\SimulationResults_Tax';
ext= '.xlsx';
filename2 = [dir '_' num2str(max_tax) ext]
g=[t_pg,fulltax,vmt_opt, opt_rev, s_rev, statetax];
xlswrite(filename2,g,1)
end



dir='D:\Projects\HeterogeneousExternalities\Output\SimulationWelfare';
ext= '.xlsx';
filename2 = [dir  ext]
z=[tlist, welfare];
xlswrite(filename2,z,1)
